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ABSTRACT 

We present an upgraded version (denoted as version 2.0) of the program 
HELAC-Onia for the automated computation of heavy-quarkonium helicity 
amplitudes within non-relativistic QCD framework. The new code has been 
designed to include many new and useful features for practical phenomenolog¬ 
ical simulations. It is designed for job submissions under cluster enviroment 
for parallel computations via Python scripts. We have interfaced HELAC- 
Onia to the parton shower Monte Carlo programs Pythia 8 and QEDPS 
to take into account the parton-shower effects. Moreover, the decay mod¬ 
ule guarantees that the program can perform the spin-entangled (cascade- 
)decay of heavy quarkonium after its generation. We have also implemented 
a reweighting method to automatically estimate the uncertainties from renor¬ 
malization and/or factorization scales as well as parton-distribution functions 
to weighted or unweighted events. A futher update is the possiblity to gener¬ 
ate one-dimensional or two-dimensional plots encoded in the analysis hies on 
the hy. Some dedicated examples are given at the end of the writeup. 



PROGRAM SUMMARY 


Program title: 

HELAC-Onia 2.0. 

Program obtainable from: 

http: //helac-phegas.web.cern.ch/helac-phegas 

Licensing provisions: none 

Operating system under which the program has been tested: 

UNIX-like platform. 

Programming language: 

Python, Fortran 77, Fortran 90, C-|--|- 

Keywords: 

heavy qnarkoninm, NRQCD, Monte Carlo simnlation 
Nature of physical problem: 

Heavy qnarkoninm prodnction processes provide an important way to investigate QCD in 
its poorly known non-pertnrbative regime. Its prodnction mechanism has been attracted 
extensive interests from the high-energy physics commnnity in decades. The qnalitative 
and qnantitative description of heavy-qnarkoninm production requires complex perturba¬ 
tive computations for high-multiplicity processes in the framework of the well established 
non-relativistic effecitive theory, NRQCD, and reliable Monte Carlo simulations to repro¬ 
duce the collider enviroment. 

Method of solution: 

Based on a recursion relation, the program is able to calculate the helicity ampltiudes 
of the high-multiplicity heavy-qurkonium-production processes. Several modules are also 
designed for dedicated simulations: 1) The code has been interfaced with the parton 
shower Monte Carlo programs; 2) A decay module to let heavy quarkonia decay with 
correct spin-correlations has been implemented; 3) The code estimates the theoretical 
uncertainties and analyzes the generated events on the fly; 4) The code is compilant with 
multi-threading/multi-core usage or cluster processors. 

CPC classification code: 

4.4 Feynman Diagrams, 11.1 General, High Energy Physics and Computing, 11.2 Phase 
Space and Event Simulation, 11.5 Quantum Chromodynamics. Lattice Gauge Theory 
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Typical running time: 

It depends on the process to be calculated and the required accuracy. 
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1 Introduction 


Since the breakthrough discovery of the Higgs boson at the LHC, much hope has been 
put on searching beyond Standard Model (BSM) particles in the next runs of LHC. 
However, the studies of quantum chromodynamics (QCD) are always playing a crucial 
role in the LHC objectives mainly because QCD is still a poorly-known theory especially 
in the non-perturbative region due to its color conhnement and because it is crucial to 
understand QCD background at the hadron colliders. For example, at 14 TeV LHC, each 
bunch crossing will generate around 50 pile-ups. Such an effect is mainly governed by soft 
interactions. Hence, all aspects of QCD still deserve to be explored as fully as possible. 

For a long time, heavy-quarkonium production and decay at high-energy colliders 
was thought to provide an ideal opporturnity to study both the perturbative and non- 
perturbative aspects of QCD. Besides, it also shows a rich physics. Fundamental parame¬ 
ters such as the strong coupling constant ag m , the heavy-flavor quark mass PEIE], the 
elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix |1], the Yukawa coupling [5] 
can be measured with heavy-quarkonium-production processes. Non-perturbative parton- 
distribution functions-either collinear ilE] or transverse [H] in the initial hadron-can be 
constrained from heavy-quarkonium data. B-meson-decay process —)■ J/i^ + pro¬ 
vides a golden channel to investigate CP violation. Heavy-quarkonium production is 
also useful in probing the multiple-parton interactions [9l [10]. Other applications in¬ 
clude quarkonium in quark-gluon plasma m, cold nuclear matter effects on quarkonium 
production [T2| and even BSM searching (see e.g. Ref. [13] and references therein) etc. 

Despite of its importance, one has very limited choice of Monte Carlo tools for the 
simulation of the heavy-quarkonium-production processes on the market. From our point 
of view, this can be attributed to several longstanding puzzles in understanding its mech¬ 
anism (see e.g Refs. [HES]) inspite of the well-established effective theory non-relativistic 
QCD (NRQCD) [IS]. Both MadOnia [T7] and HELAC-Onia [T3] are such tools dedi¬ 
cated to matrix-element calculations and event generation within the NRQCD framework, 
which aim at providing general and user-friendly public tools for theorists and experimen¬ 
talists to study the quarkonium physics. Although there are many similarities in both 
tools, we wish to emphasize some main differences between MadOnia and HELAC- 
Onia. HELAC-Onia is based on recursion relations to calculate helicity amplitudes, 
while MadOnia uses the traditional Feynman diagrams. Moreover, HELAC-Onia is de¬ 
signed to deal with processes containing one or more heavy quarkonium up to P-wave Fock 
states, while the number of states of heavy quarkonia is restricted to one in MadOnia. 

The aim of this writeup is to introduce a 2.0 version of HELAC-Onia, where many new 
and useful features are included, which are motivated by the practical phenomenological 
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simulations and the user experience, e.g. interfacing with parton shower Monte Carlo 
programs. In section we will describe the methodology, the related algorithms and 
the new features in HELAC-Onia 2.0. Then, we will show how to use the program in 
sectionSeveral examples are given in section]^ In sectionwe draw our conclusions. 
Some useful information are given in the appendices. The program structure is sketched 
in appendix]^ A summary of the new particle symbols in HELAC-Onia 2.0 is tabulated 
in appendix 1^ and a few of useful new parameters are introduced in appendix Finally, 
the addon codes in HELAC-Onia 2.0 are introduced in appendix [Dj 


2 Methodology, algorithm and new features 

2.1 Heavy-quarkonium-amplitude computation with recursion 
relations 

As it was introduced in our previous document [18], HELAC-Onia is based on the off-shell 
recursion relation [T9|. HELAC-Onia is based on a public package HELAC [201 1^ |22] . 
which is based on the Dyson-Schwinger equations [221 IMl [2S] to calculate the helicity 
amplitude in the SM at parton level. In this section, we will hrst briefly recall how to 
calculate a helicity amplitude for a general process with n external legs. We denote the 
momenta of these external legs as pi,p 2 ,... ,Pn, and their quantum numbers, such as 
color and helicity, are symbolized as ai,a 2 , ■ ■ ■ ,an- Any k external legs can form an off- 
shell current as ... ,Pi,,}] • • • > o^ifeD-^e can assign each current J' a number I, 

which is called “level”. It is dehned as the number of external legs involved in the current 
y, i.e. the “level” of y({pu,... { 0 ;*^,..., is k. The construction of the higher 

“level” currents is from the lower “level” ones in a recursion relation, where the starting 
point of the recursion relation is external legs and its end point is to obtain the “level” 
n current. The advantages by working in this way is that one is able to avoid computing 
identical subgraphs contributing to different Feynman diagrams more than once. The 
summation of all subgraphs contributing to a specihc current reduces the total number 
of objects that should be used in the next recursion procedure. 

In HELAC-Onia, we calculate the heavy-quarkonium amplitude in the framework of 
NRQCD factorization. With this formalism, the cross section for a heavy quarkonium 
Q production can be factorized into the perturbative short-distance components and the 
non-perturbative long-distance matrix elements (LDMEs). For example, at a proton- 
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( 1 ) 


proton collider, the cross section can be written as 

(t{pP ^ Q + X) = j dxidx 2 /i/p(xi)/j/p(x 2 )(j(ij ^ QQ[n] + X)(C>f), 

i,j,n 


where and are the parton distribution functions (PDFs),d(ij —)• QQ[n] +X) is the 
short distance coefficient to produce a heavy quark pair QQ in the specihc quantum state 
n. Following the usual notation, the Fock states n can be written in the spectroscopic form 
n = where S', L, J identify the spin, orbital momentum, total angular momentum 

states respectivley, and c = 1,8 means that the intermediate state QQ can be in a 
color-singlet or a color-octet state. The LDMEs are denoted as (O^). Their physical 
interpretation is a probabilitjj^ for a heavy quark pair in the Fock state n to evolve 
into a quarkonium. The power counting rules in NRQCD yield to the fact that for any 
quarkonium, there should be only a limited number of Fock states contributing to a specihc 
order of v, where v is the relative velocity of the heavy quark pair. The projection method 
is used to project a heavy-havor quark pair onto a specihc Fock state. The color-singlet 
projector is while the color-octet projector is where A“ is the Gell-Mann matrix 

and it will be projected further onto a color-how basis [261 EH EH]. The spin projectors |29| 
are in the form of 


-n(p2, A2)F5- 


+ 2E 


u{pi,\i), 


( 2 ) 


where rriQ is the mass of the heavy quark, pi,p 2 and Ai, A 2 are the momenta and helicity 
of the heavy quarks respectively. The total momentum of the heavy quark pair is = 
Pi + P 2 E = F 5 is 75 for the spin singlet state 5 = 0, and it is for the 

spin triplet state 5 = 1 with A^ = ±, 0. In order to construct the “level” 1 current for 
the heavy quarkonium, we cut the fermion chain at the place of -|- 2E in the projector 
shown in Eq.(2). Then, the new “level” I = 1 current for Q as ^u{P,X'){p + mg) and 
for Q as 


■ g72mQ ~ HELAC-Onia is also designed to be able to handle 

P-wave states. In HELAC-Onia, we introduced numerical stable P-wave currents, which 
have already been discussed in Ref. 


2.2 Angular distributions of heavy-quarkonium decays 

Angular distributions of heavy-quarkonium decays have attracted a lot of attention in 
the past few years. Theorists are interested in the polarization observables of heavy- 
quarkonium production because they might provide a “smoking gun” to discriminate 

Rigorously speaking, the LDMEs are not physical, while only the (differential) cross section is. They 
are much like the PDFs and fragmentation functions (FFs). 
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the various heavy-quarkonium-production mechanisms. The understanding of heavy- 
quarkonium polarization is also crucial for the simulations in the experimental analy¬ 
ses, e.g the detector acceptance for lepton pairs from the decay of = 1 heavy 
quarkonium strongly depends on its polarization or on the angular distribution of its 
decay products. The crue implementation of the polarization in the simulation of heavy- 
quarkonium production usually leads to one of the largest systematic uncertainties on 
the measurements. However, the current available Monte Carlo programs like Pythia are 
usually very limited in the available decay processes of heavy quarkonium and/or assume 
the unpolarized pattern. Hence, it is our motivation to implement some frequently used 
decay processes of heavy quarkonium with the polarization pattern. 

For the simple decay processes, like —)■ we only have to follow the po¬ 

larization of the mother particle and the implementation of the angular distribution 
in each spin eigenstate is straightforward. However, for a general decay process, e.g. 
Xc,j J/'tp -|- 7 , J = 1,2, the algorithms for generating the angular distributions of the 
decay products is following: 

1. Considering the helicity amplitude for the decay process is M(x), where x is the set 
of variables to characterize the kinematics of the decay process. 

2. The maximal weight of |M(x)p is H dna^ . 

3. Randomly generate a phase space point x. 

4. Uniformly generate a random number r G [0, 1]. If |M(x)p > r x hFmax, the event 
corresponding to x is retained. Otherwise, go to the former step. 

All of the available hard-coded decay processes in HELAC-Onia 2.0 can be found in de¬ 
cay/decay _list.txt. HELAC-Onia 2.0 also supports the cascade decays. For instance, 
a generated Xc,j meson can be decayed into & J/iIj and a photon 7 hrst and then the decay 
product J/^lJ can be further decayed into a lepton pair. However, the hrst decay process 
may depend on the true masses of Xc,j and J/t/. The corresponding input values can be 
changed by the user in input/decay_param_user.inp. 

Such a module is hexible and can be extended to other decay processes. Especially, 
we are planning to include EvtGen [HUl IHT] for B-physics studies in HELAC-Onia. 

2.3 Interface to parton shower Monte Carlo programs 

Parton shower Monte Carlo programs are widely used in numerical simulations for the 
collider enviroment. 
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2.3.1 HELAC-Onia+Pythia8 


Pythia is a general purpose Monte Carlo program. It provides the QCD and QED par- 
ton shower as well as the hadronization. Experiments performed on high-energy colliders 
rely heavily on it. Hence, the interface between HELAC-Onia and Pythia would surely 
extend the applications of the program. HELAC-Onia 2.0 has indeed been sucess- 
fully interfaced with Pythia^ which is written in C-I--I-. Its usage in HELAC-Onia 
requires the user to pre-install HepMC [M] and PythiaS [22] • Inheriting to its processor 
PythiaO [3S|, PythiaS provides an interface to the external hard matrix element/event 
generators via Les Houches Event hies [36| according to Les Houches accord format 071. 
In HELAC-Onia 2.0, several hies in C-|--|- are written to use the generated Les Houches 
hies and to shower and to hadronize the unweighted events with PythiaS on the hy. 
The default output is HepMC event hie after passing through PythiaS. However, such 
a format is usually inefficient to store events since it might result in a huge HepMC hie 
from a relative large Les Houches hie (say one million events). Two alternative options 
are provided. One is to output TopDrawer format plots with Hbook]^ However, such 
option requires the user to dehne all of the observables and the histograms in Fortran 
90 before calling Pythia. Useful analysis tools, like Fast Jet [3S] (or FJCore) and HEP- 
TopTagger [391 Sni E], can also be linked to hll the histograms. Some examples are 
given in analysis/P YTHIA8. Another option is using the software Root, which how¬ 
ever requires the user to pre-install Root. Events after showering and hadronization will 
be hlled into Root tree, and a pre-dehned C-I--I- Root script is necessary. We also give 
some examples in the subdirectory analysis/PYTHIA8. We will described its detailed 
usage in section 

In order to read the events record in a HepMC hie, we also provide some useful tools 
for converting it to a TopDrawer hie or a Root tree hie in HELAC-Onia 2.0. Their 
exectuable hies are HepMC2Plot and HepMC2Root in the directory bin. 

In principle, such a methodology can be applied to the interfaces to other parton shower 
Monte Carlo programs as long as it can recognize the Les Houches event hies. Although 
the current version of HELAC-Onia is still not interfaceable with other general-purpose 
parton shower Monte Carlo programs automatically, it is in our to-do list to write the 
similar interfaces for Pythia6 [33].Herwig6 [121112] and Herwig-|— f- [HI 112] . 

2.3.2 HELAC-ONiA-hQEDPS 

QEDPS is a program for the photon showering from the initial in electron-positron 

^Currently, it only works for PythiaS.I [22] but not for Pythia8.2 [55] . 

^We used a simplified version of HboOK written by M. Mangano. 



collisions [ig im sn sg, which is based on the fact that the strnctnre fnnction of an 
electron Df,±{x,Q^) shonld obey the Altarelli-Parisi eqnation 


dDe±{x,Q^) a , 2 ^ 


(3) 


where x is the longitndinal fraction, is the virtnality and P+{x) is the Altarelli-Parisi 
splitting fnnction. In the leading logarithm approximation, one can solve it to be 


a 


r-Q^ 


ds. 


-l-e 


dy 


D,4x,Q^) = U{Q\Qi)D,±{x,Qi) + — -U{Q\s) -d^P(y)D,U-^s ], (4) 


271 


>Q 
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y 


X 


where ~ ml and If is the Sndakov factor 


n(Q^Qo) = exp / —f I 

\ 27rjQ2 s Jo 


dxP{x) 


(5) 


As it is well known, the Sndakov factor has the interpretation of the probability of the 
electron evolving from scale Qq to withont emitting any hard photon. We review the 
algorithm of the photon shower in QEDPS [IHl SZl SSI SS]: 

1. Initially, set x = 1, where x is the fraction of the light-cone momentnm of the 
electron / positron. 


2. Generate a random nnmber r. If r is smaller than Il{Q^,Ql), the evolntion stops. 
Otherwise, try to hnd the next virtnality with r = Ql). At Q^, a branching 

—>■ e^(Q ^)7 happens. 


3. According to the probability of the splitting fnnction P{y) between 0 and 1 — e, try 
to determine the valne of y. Replace the original x to be xy. Iterate the step 2 nntil 
the evolntion stops. 


Hence, QEDPS provides a leading logarithm approximation to the initial state radition 
in electron-positron collisions. 

HELAC-Onia 2.0 provides an interface to QEDPS when calcnlating an eletron-positron 
annihilation process. Thanks to the release of the sonrce code, QEDPS is self con¬ 
tained and the nser does not need to install it by himself/herself. An inpnt parame¬ 
ter emep_ISR_shower is provided in input/user.inp or input/default.inp to determine 
whether nsing QEDPS to perform a photon shower (see its description in appendix]^. If 
one wants to generate TopDrawer, Gnuplot or Root plot hies, one shonld also edit the 
Eortran snbrontine plot_f ill.QEDPS before compilation. We have performed an appli¬ 
cation of HELAC-Onia-|-QEDPS to J/Jj inclnsive prodnction at B factories in Ref. [50] . 
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2.4 Estimating the scale and PDF uncertainties 

Varying the renormalization scale fiR and the factorization scale /xj? is often thought to be 
a standard way to estimate the theoretical uncertainty in perturbative computations due 
to the missing higher-order contributions. Although such an argument can be applied to 
the scattering or decay processes in general as long as its validation of the perturbative 
description, there are indeed several cases which we already encountered where it is not 
applicable, such as where one encounters large coefficients correction (from large loga¬ 
rithms, large tt^, or large color factors) or new channels (e.g. new initial states, new phase 
space region, or new fragmentation topology). Unfortunately, the later case frequently 
happens in heavy-quarkonium-production processes. Because new fragmentation topolo¬ 
gies appear only at higher orders in perturbative calculations, it is usually necessary to 
consider the contributions from the higher-multiplicity processes accompanying with more 
partons. Some examples indeed already can be seen in the single-quarkonium |511152] and 
double-quarkonium [53] production processes. Nevertheless, after taking into account all 
of the important topologies, the scale dependence is sufficiently reasonable to estimate the 
higher-order corrections. At tree level, the renormalization-scale dependence is only in 
the renormalization running of ctsif^R), while the factorization-scale dependence is in the 
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution |5H (55] [56] of the PDFs. 
Hence, it is straightfoward that the estimation of the scale uncertainties are irrevelant 
to the most time-consuming matrix-element calculations as long as one knows the initial 
states and the perturbative orders. From the technical point of view, such estimation can 
be zero GPU costj^ In calculating a physical observable or hlling a histogram, one just 
multiplies a weight 




to the central value in each phase space point or each event, where fi is the PDF, x* is 
the Bjorken fraction, (central) renormalization and factorization 

scales, b is the power of a* in the squared amplitude. Such reweighting procedure has 
been widely used in other programs, such as MadGraph [58l (59] EO] . 

Another important source of theoretical uncertainty that can be obtained from the 
reweighting method is the PDF uncertainty, which does not reflect the uncertainty in the 
hard matrix element but rather the uncertainty in the extraction of the PDF. It is known 
that the PDF uncertainty stems from at least three sources: the uncertainties in the in¬ 
put (experimental) data, the accuracy of the perturbative calculation, and the method to 

^In a general case, the explicit renormalization and factorization scales dependence is also known at 
next-to-leading order (NLO) [57| . 
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extract the PDF. Most of the modern PDF sets provide a way to estimate the impact of 
such uncertainty to the theoretical calculations. For example, the global-fit PDF MSTW 
2008 NLO ini] provides 40 error PDFs to quantify its uncertainty. Instead of reevaluat¬ 
ing the matrix element with new PDF, one is able to evaluate the PDF uncertainty by 
multiplying a weight 


Wpdf(/(,/2;/i,/2) 


fl{f^F,Xi)f2{fiF,X2)' 


(7) 


where is an error PDF and fi is the central PDF. Such a procedure is exact in a 
parton-level calculation. However, it should be understood as an approximation when 
incorporating with a parton shower Monte Carlo program, since the backward evolution 
of the initial state partons in Monte Carlo indeed contains an implicit dependence on the 
chosen PDF, i.e. the central one in a Les Houches event hie. 

Although such a procedure is more trivial at tree level than that at NLO [57], for 
completeness, we would like to emphasize that HELAC-Onia 2.0 provides a functionality 
to estimate scale and PDF uncertainties with such a reweighting method. However, 
because of the recursion relations, it is non-trivial to separate different coupling orders 
without degrading its speed advantage. Hence, we want to remind the user that it would 
be wrong in evaluating renormalization scale dependence if the amplitudes in different 
orders will contribute to a partonic level process. 


2.5 A summary of new features 

We give a summary before closing this section. HELAC-Onia 2.0 has been improved 
much compared to the hrst released version of HELAC-Onia [18]. The main changes 
include: 

1. Two completely independent generators based on PHEGAS [62] and VEGAS [63] 
are implemented. Both of them can generate unweighted events for 2 —)■ u processes 
when u > 2 at pp,pp and e“e+ collisions. For 2 —)■ 1 processes at hadron colliders, 
only VEGAS is available. 

2. Additional internal PDFs are included. The program can also be interfaced with 
LHAPDF [61] • 

3. Analysis plots are done on the fly. Differential distributions can be plotted at the 
end of the computation phase. 

4. The laboratory frame is not restricted to the center-of-mass frame of the initial 
collision anymore. The fixed-target mode was also added. 
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5. An interface from QEDPS to HELAC-Onia is done. One can include the QED 
photon showering effects from initial beams. 

6 . An interface from Pythia 8 to HELAC-Onia is done. It is able to use the major 
functionality in Pythia 8 . 

7. Reweighting method is used to estimate the renormalization/factorization scale and 
PDF uncertainties. 

8 . Several spin-entangled decay processes are implemented to take into account the 
polarization effects. 

9. A user-friendly Python script is provided for user to use the program. It also allows 
us to calculate the cross sections of several subprocesses with multiply CPUs, such 
as on a multicore computer or on a cluster. 

For item 1, in the previous version, the unweighted events can only be generated by 
PHEGAS [62]. 2 —)• 1 processes are not handled in this case. This improvement allows us 
to lift several restrictions. Concerning item 2, beforehand, only CTEQ6 [65] was available. 
It paves the way to application in nucleus collisions and to estimate PDF uncertainties. 
With the help of the improvement presented in item 4, we are able to apply HELAC- 
Onia to more experiments like hxed-target experiments [HS] • Item 5 allows us to consider 
initial radiation effects at collision, which might not be small in several important 
processes,such as —)■ J/ip + gg [50], while the improvement of item 6 extends the 

usage of HELAC-Onia to one of the most widely used multipurpose event generator 
Pythia 8 . For item 7, it will be very useful to estimate the renormalization/factorization 
scale and PDF uncertainties without the extra cost of recalculating cross sections. The 
improvement of item 8 is quite useful for practical simulations, and that of item 9 improves 
the user experience on using the program. 

3 How to use the program 

In this section, we will hrst give a brief introduction on how to perform a phenomenological 
analysis in a basic way. If one is only interested in using the program, one can follow the 
instruction in this section and ignore the remaining context. 

3.1 Standalone 

We hrst introduce how to use HELAC-Onia 2.0 in a standalone way. Before running the 
code, one should specify the conhgurations via the conhguration hie ho_configuration.txt 
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in the input subdirectory, which is described in section For example, if one wants to 
use LHAPDF [M] , one should assign the correct path to the parameter lhapdf in the hie 
ho_configuration.txt. Useful comments are also given in the conhguration hie. If one 
wants to output plot(s) on the hy, one should also edit the user plot hie plot_user.f90 in 
the subsubdirectory analysis/user. One can follow some example hies to write his/her 
own plot hie. After the above preparations, one can set the conhguration and make the 
hies via the command line 

I > ./config 

This procedure should only be done once. Afterward, the program is ready for running. 

In HELAC-Onia 2.0, one can use two modes to perform a computation of a cross 
section. We still keep the initial way to run the program directly via exectuable hie 
Helac-Onia. To use this way, one should follow the following lines: 

1. Specify input parameters in input/user.inp following the format in input/default.inp. 
Some examples are given in input/demo. 

2. Provide the process information in input/process.inp as well as the decay infor¬ 
mation in input/decay user.inp. 

3. If one wants to dehne one’s own dynamical renormalization and/or factorization 
scale, one should edit it in src/setscale.f90 before compiling. Four default scales 
are dehned, i.e. the hxed scale, the transverse mass rriT^i = \Jrn\ + of the hrst 

hnal state, Po = y^(E/efinal states rrifY + Pl^ 5 E/efinal states + 

4. Run the program with the command line 

I > . /Helac-Dnia 

or 

I > ./bin/Helac-Onia 

The hnal results will be collected in the output directory. 

A second way to run HELAC-Onia 2.0 is by using the Python scripts, which is more 
user-friendly and hence strongly recommended. It provides the opporturnity to avoid 
mixing the working directory and the HELAC- Onia directory. If one wants to use the 
majority of the new features in the program, one has to run the program with the Python 
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scripts. The basic way of using it after the above hrst three items is to run the program 
with the script 

I > /PATH/TO/HELAC-Onia/ho_cluster 

where /PATH/TO/HELAC-Onia is a path to the HELAC-Onia directory relative to 
your working directory. Then one will see a prompt starting with “HO>”. There are 
two phases to compute a cross section for a process, i.e. generate the process and run the 
program. In the hrst phase, one should dehne one process or several subprocesses. For 
example, if one wants to calculate J/'ip pair production, the syntax should be 

I H0> generate g g > cc~(3Sll) cc~(3Sll) 

where the symbol g represents the initial gluon and cc ~ (SS'll) means a pair of charm 
(anti-)quark c and c in fi conhguration, which can be found in appendix]^ One can also 
calculate the cross sections for several subprocesses simultaneously. For instance, if one 
wants to go beyond the leading-order computation of J/V’ pair production, he/she can 
type the following command lines 

H0> generate g g > cc~(3Sll) cc~(3Sll) 

H0> generate g g > cc~(3Sll) cc~(3Sll) g 

H0> generate u g > cc~(3Sll) cc~(3Sll) u 

H0> generate g u > cc~(3Sll) cc~(3Sll) u 

It will use multiple cores to calculate the cross sections on the cluster or on a multicore 
machine. One can also run the addon processes, where the available addon processes and 
the corresponding numbers are listed in addon/addon_process.txt. For instance, if one 
wants to calculate the double parton scattering (DPS) for J/'ijj pair production, one can 
generate the process via 

I H0> generate addon 1 

where the keyword addon should be specihed after generate and the number for this 
process is 1 as be seen in addon/addon_process.txt. Before launching the jobs for 
numerical calculations, one can also change the parameters in input/user.inp via the 
interactive command syntax 

I H0> set <parameter_name> = <value> 
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One example is to take the maximum Monte Carlo integration number to be 10000. Then 
one just simply types 

I H0> set nmc = 10000 

If one wants to take VEGAS as the Monte Carlo integration program, please uses 
I H0> set gener = 3 

Another useful feature is to dehne the decay process(es) via Python scripts. The syntax 

is 

I H0> decay <process> @ <branching ratio> 

The command lines 

H0> decay cc~(3Sll) > ni+ m- 0 0.06d0 
H0> decay w+ > m+ vm 0 IdO 

means that the fi charmonium in the hnal states will decay to a muon pair with the 
branching ratio 6% and the W~^ boson will perform leptonic decay to a muon and a 
neutrino with 100% probability. After all, one just submits the job via command 

I H0> launch 

and waits for the hnal results to be collected in a new created subsubdirectory PROC_HO J/results 
in the working directory, where i is a number to be assigned uniquely. 

In the new way, we take a similar fashion of the widely used program MadGraph5_aMC@NLO, 
and we hope that it would become a standard in the future, or at least it will be much 
easier for the users who are already familiar with MadGraph5_aMC@NLO to learn to 
use this program. 

3.2 HELAC-Onia+Pythia8 

Let us start to consider the case of using HELAC-Onia+PythiaS]^ One has to gener¬ 
ate the Les Houches hie for unweighted events in PROC_HO_i/PO_calc_j/output hrst 
before calling Pythia. This can be achieved by setting the hag unwgt to be T before 
launching the program. 

®It is more or less trivial to use HELAC-OniA-|-QEDPS by setting the flag emep_ISR_shower to be 
1 in input/user.inp as described in appendix 
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The program will be able to call PythiaS if one sets correct path for PythiaS and 
HepMC in the conhguration hie ho_configuration.txt and specihes the shower param¬ 
eter^ in the hie shower _card_user.inp before using the ho .cluster script. The corre¬ 
sponding parameter setup for calling PythiaS in user.inp is 

I H0> set parton.shower = 1 

All of the above shower related setup can be done before or after generating Les 
Houches event hies. In the former case, the PythiaS will be called with the command 
launch directly, while in the later case we also provide a new command shower with the 
syntax 

I H0> shower <working path> 

where the working path is usually PROC_HO J. The hnal output from PythiaS will be 
collected in the directory PROC_HO J/PO_calc_j/shower/HO_PYTHIA8_k, where 
i,j and k are integers starting from 0. 

We will present some technical details for using PythiaS in HELAC-Onia 2.0 directly, 
since the user might have encountering problems in compiling the PythiaS related code 
with wrong setup. It will also be useful for extensive usage of PythiaS in HELAC-Onia. 

After running the program with launch or shower commands with the correct shower- 
related setup, a directory HO_PYTHIA8_k will be created in PROC_HO_i/PO_calc_j/shower. 
If the program compiles successfully, an exectuable hie Pythia8.exe will be generated. 

One can also change the PythiaS setup via the hie Pythia8_lhe.cmnd in the same 
directory. A direct using PythiaS is possible by simply typing 

I > . / PythiaS.exe 


If there does not exist the exectuable hie Pythia8.exe, it means there is a problem in the 
compilation. Some useful information can be found in shower.log to solve the problems. 


^Especially, one should be aware of the parameter ANALYSE in shower_card_user.inp. It determines 

If ANALYSE is empty, it will output a HepMC event 


the output mode as described in section 2.3.1 


file from Pythia. One can also take ANALYSE to be a PortraN90 or C-|—h plot file in the other 
two output modes in section |2.3.1[ For example, if ANALYSE = plot_py8_pp_t j and there exist a file 
plot_py8_pp_tj .f90 in analysis/PYTHIAS in the HELAC-Onia root directory, it will output a 
TopDraWER file from PythiaS. If the extension of the file is . cc in analysis/PYTHIAS, it will 


output a Root file instead. 
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4 Examples 

4.1 NNLO* level J/'ip and 'ip{2S) hadroproduction 

-0 and T production at hadron colliders have challenged our understanding of heavy- 
quarkonium mechanism for decades in EH Eg. Since then, the heavy-quarkonium- 
production data have been removed in the global £t of extracting PDF. Large QCD 
corrections were found in heavy-quarkonium production due to new p'^-enhanced frag- 
mentating Feynman diagrams at higher orders lai. Hence, it was suggested to look at 
how partial next-to-next-to-leading order (NNLO) QCD correction impacts the differen¬ 
tial cross sections of V’ and T production [52], which was called NNLO*. Besides of the 
complete NLO result [SI], it requires to calculate the 0{a^g) tree-level 2 —)■ 4 process 
pp — )■ -^-l-S-jets. MadOnia mi was able to perform a hrst numerical computation with 
such complex. It is a good process to show the robustness of HELAC-Onia and to 
compare MadOnia. 

The calculation of 0(0;^) process pp —)■ "0-1-3-jets in the color-singlet mechanism (CSM) 
consists the following 13 independent subprocesses: 


H0> 

generate 

g 

g 

> 

cc" (3S11) 

g 

g 

g 

H0> 

generate 

g 

g 

> 

cc" (3S11) 

u 

u~ 

g 

H0> 

generate 

u 

g 

> 

cc" (3S11) 

g 

g 

u 

H0> 

generate 

g 

u 

> 

cc" (3S11) 

g 

g 

u 

H0> 

generate 

u 

g 

> 

cc" (3S11) 

u 

u~ 

u 

H0> 

generate 

g 

u 

> 

cc" (3S11) 

u 

u~ 

u 

H0> 

generate 

u 

g 

> 

cc" (3S11) 

d 

d~ 

u 

H0> 

generate 

g 

u 

> 

cc" (3S11) 

d 

d~ 

u 

H0> 

generate 

u 


> 

cc~ (3S11) 

g 

: g 

g 

H0> 

generate 

u 

u*" 

> 

cc~ (3S11) 

u 

. u 

g 

H0> 

generate 

u 


> 

cc~ (3S11) 

d 

. d 

g 

H0> 

generate 

u 

u 

> 

cc~ (3S11) 

u 

u 

g 

H0> 

generate 

u 

d 

> 

cc~(3Sll) 

u 

d 

g 


In the above command lines, we only include up and down quarks. The other quark- 
initial-state contribution can be included by 

H0> set quarksumQ = T 
H0> set iqnum = 3 

The truth of the flag quarksumQ makes sure we will include the initial (anti)quark PDFs 
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and iqnuni=3 means we are working in 3-light-qnark-flavor scheme. Moreover, one can 
type 

I H0> set combine_f actors = 1. 3. 1. 1. 1. 1. 2. 2. 1. 1. 2. 1. 1. 

to explicitly mnltiply a combination factor for each snbprocess. For example, in the snb- 
process gg —>■ cc^]+uug, we take a factor 3 to acconnt for gg —>■ cc^]+qqg with q = u,d, s, 
becanse they share the same matrix element and PDF. The detailed correspondances are 
shown in Tabj^ In order to avoid infrared divergence in NNLO* calculations, a special 
cutoff [52] should be applied to the invariant mass squared of any massless parton 
pair, i.e. {pi + PjY > In this case, we set 

I H0> set minmqqp = 3d0 

for any hnal state massless parton pair a /{pi + PjY ^ 3.0, and set 
I H0> set minmqbeam = 3d0 

for one hnal state massless parton and a initial state parton — PiY — 3-0. 

We have compared the HELAC-Onia result with the MadOnia result, and found they 
were in perfect agreement. Because HELAC-Onia is based on the recursion relations, 
HELAC-Onia is faster than MadOnia in the computations. It is much easier for us to 
extend the MadOnia result to a wider pt range. In Ref. coi. ATLAS Collaboration 
already used our result to compare their measurment for 'ip{2S) up to 100 GeV. Here, 
we present the px distributions of J/ip (FigQ in the LHCb acceptance at 13 TeV. We 
should remind the reader the following points. The color-singlet LDME is estimated in 
potential model EH. The corresponding radial wave function at the orgin was derived 
in the QCD-motivated Buchmuller-Tye potential [72||^ We used CTEQ6M [65] as our 
PDF set and hxed = 2mc- The error bands in Figj^ represent the renormalization 

and factorization scale uncertainties +V t_ < < 2^/{2/mPp^/^^ and the 

uncertainty in charm quark mass nic = 1.5 ±0.1 GeV. In Figj^ we do not include the 
NLO contribution, although it should be contained for a real NNLO* prediction. 

4.2 pp ^ J/t/) + J/t/) + cc 

As noticed in Ref. [TH], in double J/ip production, pp —>■ J/ip + J/ip + cc shares the 
leading-pT contribution from charm quark fragmentation diagrams though it is of 0{al) 

^We will use the same color-singlet LDME in the following computations. 
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Syntax 

Subprocess 

g g > cc~(3S11) g g g 

gg cc[6] + ggg 

g g > cc~(3Sll) u u~ g 

gg -)■ cc[B] + qqg 
with q = u,d,s 

u g > cc~(3Sll) g g u 

qg cc[B] + ggq 
with q = u,u, d, d, s, s 

g u > cc~(3Sll) g g u 

gq cc[B] + ggq 
with q = u,u, d, d, s, s 

u g > cc~(3Sll) u u~ u 

qg —)■ cc[B] + qqq 
with q = u,u, d, d, s, s 

g u > cc~(3Sll) u u~ u 

gq —)■ cc[B] + qqq 
with q = u,u, d, d, s, s 

u g > cc~(3Sll) d d~ u 

qg —)■ cc[B] + q'q'q 
with q = u,u, d, d, s, s 
and q' = u,d,s 

and q, q' not in the same flavor 

g u > cc~(3Sll) d d~ u 

gq —)■ cc[B] + q'q'q 
with q = u,u, d, d, s, s 
and q' = u,d,s 

and q, q' not in the same flavor 

u u~ > cc~(3Sll) g g g 

qq cc[B] + ggg 
with q = u,u, d, d, s, s 

u u~ > cc~(3Sll) u u~ g 

qq —)■ cc[B] + qqg 
with q = u,u, d, d, s, s 

u u~ > cc~(3Sll) d d~ g 

qq —)■ cc[B] + q’q'g 
with q = u,u, d, d, s, s 
and q' = u,d,s 

and q, q' not in the same flavor 

u u > cc~(3Sll) u u g 

qq —)■ cc[B] + qqg 
with q = u,u, d, d, s, s 

u d > cc~(3Sll) u d g 

qq' —)■ cc[B] + qq'g 
with g, q' = u, it, d, d, s, s 
and q, q' not in the same flavor 


Table 1: Subprocesses are calculated with^^ch generation for pp —)■ '^+3-jets in CSM 
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Figure 1: The transverse momentum distributions of J/^|J from pp —)■ J/ijj + jjj in the 
LHCb acceptance at 13 TeV. 

suppressed compared to pp —)■ + Hence, it is necessary to quantify its maginitude 

by an explict calculation. For the hrst time, we performed such a complex calculation 
with the help of HELAC-Onia 2.0 in Ref. [lO], which involves more than 2000 Feynman 
diagrams. It is a hrst ( and till now the only ) 2 —)■ 4 process with at least two quarkonia 
to be calculated. We take this example to show the uniqueness and the rubostness of 
HELAC-Onia to perform perturbative computations of more than one quarkonium pro¬ 
duction processes. Because the luminosity of the quark-antiquark initial states is usually 
much smaller than that of the gluon-gluon initial state at the high-energy colliders, we 
only include the gluon-gluon initial state here. One can use the following command line 
to generate the process 

I H0> generate g g > cc~(3Sll) cc~(3Sll) c c~ 
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For illustration, we are working in the CMS acceptance |73|: 


p^J^ > 6.5 GeV 

if 

|//^| < 1.2, 


> 6.5 ^ 4.5 GeV 

if 

1.2 < < 1.43, 


p-!/^ > 4.5 GeV 

if 

1.43 < < 2.2, 

(8) 


where in the rapidity interval 1.2 < < 1.43, the transverse momentum cutoff 

scales linearly with its absolute rapidity We used the exact setup in Ref. HD], 

Some selected differential distributions inpp—)■ JI'll) ^ Jjil) ^ cc are shown in Figj^ The 
error bands are coming from the variations of the renormalization and factorization scales 
and the uncertainty of charm quark mass. Following the way in Ref. [10], we have taken 
into account the feeddown contribution from 'ip{2S) decay. The feeddown contribution 
enhances the (differential) cross section by a factor of 1.89. The absolute azimuthal 
difference between the J/ip pair is shown in Fig, 2a It is an observable to distuiguish 
the double-parton scattering (DPS) and the traditional single-parton scattering (SPS) 
since in the former production mechanism the two J/ip are uncorrelated. However, such 
an observable might be sensitive to the primordial smearing from the beam 


Besides the absolute azimuthal difference, the absolute rapidity difference (Fig, 2b) and 
the invariant mass distribution (Fig, 2c) are also good kinematical variables to discriminate 
DPS and SPS. Finally, various transverse momentum spectra are displayed. In Figj2d[ we 
presented distribution of the vectorial transverse momentum sum 


2e 

T/p 


(Fig,2f) is the yileds of the leading pt = (subleading px 


while Fig 

of fhe two J/ip. 


4.3 J/iIj hadroproduction with parton shower effect 

The inclusive J/ip hadroproduction is a hrst process challenging our understanding of 
the heavy-quarkonium-production mechanism. For a long time, it was known that the 
CSM can describe the total cross section of J/ip or T (e.g. see Ref. [71|) but not in 
the transverse momentum px distributions. In the recent years, most of the studies were 
focusing on the interpretation of the yields [SH E21 EH EH EH EHl EH EH EH E2] and the 
polarization [SH EH EH EH EH EH EH EH EH E2] of single-quarkonium hadroproduction 
at the large px regime. Some efforts have also been paid in the small px regime [HH |9l] . 
However, none of the consistent matching between large px and small px results exists. 
In Ref. [HH , analytical small px resummation is performed for the color-octet states only 
in NRQCD, which lacks the dominant color-singlet contribution and the matching to the 
hxed-order results. Alternatively, one can perform a resummation with the parton shower 
(PS) approach, which is formally to be restricted to the leading-log accuracy (although 
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the partial subleading-log contributions can also be taken into account). It generates the 
complete events with correct kinematics and can be applied directly on the experimental 
analysis by including the detector effect. 

In this subsection, we will give a simple example to show the importance of parton 
shower effect for J/'il) hadropro duct ion in the small pr region. Let us consider the color- 
singlet contribution only at leading order (LO) in as and in v (the relative velocity between 
the charm quark pair). Without primordial kx smearing effect from the beam and the 
multiple interactions, the LO curve in px distribution is indeed siginificantly smearing by 
PS as seen in Figj^ where we have used Pythia8.186. In the left pannel of Figj^ one 
can observe that such a smearing effect is mainly from the initial state radiation (ISR) 
while the final state radiation (FSR) only distorts the distribution slightly. On the other 
hand,LO result is good enough to describe the rapidity distribution. We should remind 
the reader that while the LO-|-P^ color-singlet contribution is expect to describe the 
small px data, the intermediate and the large px data will receive substaintial higher- 
order (or real emissions) [SHEg and color-octet contributions. A consistent treatment of 
J/'ip production in NRQCD is possible with the LO merging of matrix elements and PS 
in different jet multiplicities ES], where a pioneer work has been done in rg, produc¬ 
tion 1^3 • Such a detailed analysis is of course interesting but beyond the scope of this 
paper. We will leave it for a future work. 


4.4 Validation of decay angular distributions 


In subsection |2.2[ we have discussed the implementation of the angular distributions in 
the heavy-quarkonium decay in HELAC-Onia. Here, we will give three examples to 
validate our implementations: —>• , Xci J/ip + 'y and Xc 2 —t J/ip + 7 . 

In our first example J/'ip ^ we know the angular distribution of one lepton is 

in the form of 


da 


dcosO 


1 -I- Ag cos^ 6*, 


(9) 


where 9 is the polar angular respect to the spin quantization axis in the rest frame of J/V’ 
and \q can be expressed in the longitudinal polarized cross section ctl and the transverse 
polarized cross section ctt of J/'tp 

_ (Tt — 2(Tl 

— --w— 

(Tt + 2(Tl 

®One should turn on primordial kx in Pythia as well 


( 10 ) 
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We have compared the numerical result from HELAC-Onia and the analytical result 
Eq.Q in Figj^ The total cross section has been normalized to unity. We selected 
Xg = 1, —1, 0, 0.3, 0.5 for illustration. Perfect agreement is found. 

For Xci —t J/ip + 7 , we have the same polar angular distribution Fq.([^ with respect 
to the decay product J/ip or x [HE], while for Xc 2 —>■ J/ip + '^ the general formula is 

, ^ 1 + Ae cos^ 6 + A26> cos^ 6 . (11) 

a cos 0 

In the later case, A 20 is suppressed by the higher-order multipole amplitudes [HE]- Explic¬ 
itly, we have for Xci 


Xg = {l- 35) 


cr, 


Xci 

tot 


3^0,0 


(1 + 6)aP/;P + (1 - 35)aJ;c 


( 12 ) 


and for Xc 2 


Xg = 6 [(1 - 35o - 5i)(Jt^of - (1 - 75o + - (3 - 5o - 75i)crJ_o] /R, 

A 20 = (1 + 55o - 55i) [appp - 5{appp + /R, 

R = (1 + 55o + 35i)o'^‘f -l- 3(1 — 35o — + (5 — 75o — 95 i)(Jq. (13) 

In the above equation, we denote is the (f, j)-component of the spin density matrix of 
Xci, 2 production. In the following, we take cr/p/’X = which is valid in a CP-conserved 

process. The spin-summed cross sections can be expressed as 


a, 


Xcl _ 


tot 


= a 


Xcl 

1,1 


I ^A.ci I XCL 
^ *^0,0 *^-1,- 


cr, 


tot 


= a. 


2,2 


1 ) 


+ + CX 


0,0 


'- 1,-1 


+ cr. 


- 2 ,- 2 - 


(14) 


Parameters 5,5o!5i enter into Eq.(12) and Eq.(13). They can be determined by the nor- 
malizec^ multipole amplitudes 


5 

50 

51 


(1 -I- 2af“^a2“^) 

2 ’ 

1 -|- 2a{^‘^{\/5a2^‘^ + 203 "^) -|- + t/Sog^^) -|- 3 

10 ’ 

96 af=^(\/ 5 a 2 ^^ — 403=^) — 402^^(03^^21/503=^)-h 7 (03=^)^ 

-4b-> (15) 


where a(~^,a 2 ~^ and a'^~^ are the electric dipole (El), magnetic quadrupole (M2) and 
electric octupole (E3) amplitudes for Xcj- We take the measured values by CLEO col¬ 
laboration in Ref. [98]. The numercial values are shown in Tabj^ However, the imple¬ 
mentation of the cascade decay Xc J/ipx —t in HELAC-Onia requires its full 


®We have (af“^)^ + ( 02 ”^)^ = 1 and -I- ( 02 ”^)^ + (^ 3 ”^)^ = 1- 
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J=J 

J=j J=j 

0-2 «3 

j 

= 1 

0.998 

-0.0626 

j 

= 2 

0.996 

-0.093 0 


Table 2: The normalized multipole amplitudes of Xcj J/V’ + 7 from the CELO mea- 
surment [98]. 


knowledge of the helicity decay amplitudes in terms of multipole amplitudes. Its com¬ 
plete derivation was performed in Ref. |99| at the amplitude level for the first time. Such 
amplitude-level experssions will be served as the helicity amplitude ^(x) defined in sub- 
The validation of the implementations for Xci,2 —t J/V’ + 7 in HELAC-Onia 


section 2.2 


2.0 can be found in Figj^and Figj^ The histograms of the decay product Jip’s angular 
distributions perfectly agree with the analytical experssions. 


5 Conclusions 

We have presented a version 2.0 of HELAC-Onia with several important updates for 
the practical theoretical studies and the Monte Carlo simulations for heavy-quarkonium- 
production processes. The main improvements are 

• a completely new interface for talking between the user and the program written 
in Python scripts. It is much user-friendly and suitable to submit calculation jobs 
with multi-threading usage or on a cluster; 

• automated interfacing HELAC-Onia to the parton shower Monte Carlo event gen¬ 
erators. Two parton shower programs are sucessfully linked; One is QEDPS for the 
initial photon showering from the processes in electron-positron collisions, while the 
other one is the widely used one Pythia 8; 

• a decay module for perfoming the spin-entangled (cascade-) decay of heavy quarko- 
nium. Some dedicated decay processes are implemented such as J/il) —)■ 

Xc ^ J/i’ + 1 and the decays of top quark, W-boson, Z boson; 

• a reweighting method for estimating the uncertainties from the renormalization/- 
factorization scale and PDF in an automatic manner; 

• one-dimensional or two-dimensional histograms generation on the fly. Moreover, we 
also provide several useful analysis tools. 
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All of the above improvements are quite useful in the study of the heavy-quarkonium 
production. It also provides a flexible framework for the future developments like heavy- 
quarkonium production in heavy ion collisions Hffll or in the transverse momentum fac¬ 
torization framework [S]. 
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A Program structure 


In this section, we will describe briefly the new program structnre of HELAC-Onia 2.0 for 
the future developments. The main files contained in the program are already described 
in the README file of the tarball. The files in the program are mainly included in 
several subdirectories, which are displayed in Figj^ There are mainly ten subdirectories 
under the main directory of HELAC-Onia. Let us explain them in somewhat detail: 

• input. All of the input files that required by the program are contained in this 
subdirectory. They are: 

— user.inp: a file for user to specify his/her input parameters. 

— default.inp: a file that includes all of the default values for the input param¬ 
eters. 

— process.inp: a file for user to tell the program the process information. 

— ho_configuration.txt: a configuration file for HELAC-Onia. 

— seed.input: a seed for random number generator. 

— shower_card_user.inp (shower_card_default.inp): a user (default) card to 
use parton shower programs. 

— decay par am user.inp (decay param default, inp): a list of user-defined 
(default) parameters for using in the decay module. 

— decay_user.inp (decay_default.inp): a file to specify decay chains in this 
card. 

• output. All of the output files will be generated here. Initially, it is empty. 

• src. It contains all of the main source files of the program. They can be mainly 
divided into two parts. One part is for the matrix elements generator and the other 
part is for the phase space integration and events generation. 

1. matrix elements generation. 

— Helac_Global.f90: It is a file which contains all of the global variables. 

— Helac_Func_l.f90: In it, many helper functions and subroutines are de¬ 
fined. 

— alfas Junctions.f90: Running of as which is used in MCFM |101] . 

— Projectors.f90: It is a file in which the Clebsch-Gordan coefficients are 
defined. 


26 




— Constants.90: Several subroutines are defined for reading input param¬ 
eters. 

— SM_FeynRule_Helac.f90 ;It contains all of the Feynman rules of the 
Standard Model. 

— Feynman_Helac.f90: A useful subroutine is written in this file for recon¬ 
structing all Feynman diagrams. 

— Helac_wavef.90:It is a file to define all of external wave functions. 

— Helac pan2.f90: Definition of vertices to be used in Helac_panl.f90. 

— Helac _panl.f90 ; Off-shell currents generation by using recursion relation. 

— Helac_master.f90: It is a main file of computing helicity amplitudes. 

2. phase space integration and events generation. It is based on several adapted 
Monte Carlo integration programs. 

(a) PHEGAS: 

— Phegas.f90 :It is an extensive version of PHEGAS [62] to deal with 
quarkonium kinematics. It was rewritten in Fortran 90. 

— Phegas_Choice.f90: Some helper functions are defined here that will 
be used by Phegas.f90. 

(b) VEGAS: 

— MC_VEGAS.f90: A Fortran 90 version of VEGAS [63] . 

— Func_PSI.f90: Some helper functions of phase space integration were 
written in this file. 

— Colliders_PSIl.f90: Phase space integration with VEGAS for 2 —)■ 
n(n > 1) at hadron colliders. 

— Colliders_PSI2.f90: Phase space integration with VEGAS for 2 — 
n{n > 2) at electron-positron colliders. 

(c) MINT: 

— mint-integrator.f90: It is a Fortran 90 version of Mint [102] . 

(d) Internal Fortran 90 PDF files: 

— CTEQ6PDF.f90:GTEQ6 PDF [65] file in Fortran 90 version. 

— Structf_PDFs.f90: A file for calling PDFs. 

(e) LHAPDFfile: 

— Structf_LHAPDF.f90: A file for calling PDFs from LHAPDF [5^ . 

User should specify “lhapdf=/path/to/lhapdf-config” in input/ho_configuration.txt. 
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(f) Others: 

— Helac_ranmar.f90; A random number generation program Ranmar 
in Fortran 90. 

- MC_PARNI_Weight.f90: PARNI in Fortran 90, but it is not used. 

- MC_RAMBO.f90: Rambo |103] in Fortran 90. 

- MC_Helac_GRID.f90: A grid file. 

— Helac_unwei.f90: There are some subroutines for dealing with un¬ 
weighted events in this file. 

— ADAPT.f90: It is for optimization by using adaption procedure. 

— Phegas_Durham.f90: Durham in Fortran 90. It can only be used 
to generate phase space points for massless external particles. 

— MC_Func.f90; There are some helper functions and subroutines for 
Monte Carlo integrations. 

— Kinetic_Func.f90: Some kinematical variables are dehned in this file. 

— Cuts_Module.f90; It is a file to provide the user to impose kinemat¬ 
ical cutoff. 

— KT_Clustering.f90: kx clustering and reweight factor for MLM merg¬ 
ing [Ml 196]. 

— setscale.fOO: It provides the user to specify his/her renormalization 
and factorization scales. 

— setscale_default.f90: It is only a default setscale.f90 file for backup. 

— Helac_histo.f90: Histogram drawing hie in HELAC. 

— SinglePro.f90: It is the main hie for phase space integration and 
events generation. 

— Summation_Pro.f90: A hie for the summation mode, which is not 
used yet. 

— unweight Jhe.f90: A hie for writing out Les Houches events hies. 

— FO_plot.f90: A hie for plotting hxed-order distributions. In this case, 
unweight events generation is not necessary. 

— Main_Test.f90; The Fortran 90 main program. 

pdf. More extensive internal PDFs are located in this subdirectory. 

— pdf_list.txt: A summary of internal PDFs in HELAC-Onia. 






— make_opts,makefile_pdf: Files of makefile for the PDF related routines. A 
library libpdf.a will be generated in lib subdirectory. 

— opendata.f; A file in Fortran 77 for opening PDF data. 

— Partonx5.f: Standalone Fortran 77 Partonx function. 

— CTEQ files: They include cteq3.f,Ctq4Fn.f,Ctq5Par.f,Ctq5Pdf.f,Ctq6Pdf.f. 

— MRST hies: They include mrs98.f,mrs98ht.f,mrs981o.f,mrs99.f,mrst2001.f,jeppe02.f. 

— gsdpdf hie: They include GS09 dPDF hies |in5] . 

• shower. The subdirectory contains hies for parton shower. 

— QEDPS: It contains the hies of QEDPS for ISR photon shower form initial 

beams. 

— PYTHIA8: Pythia 8 subsubdirectory. It includes the main hies for interfac¬ 
ing HELAC-Onia to Pythia 8 for showering. 

— PYTHIA6: Pythia 6 [35] subsubdirectory. It will be used for the future 
development. 

— HERWIG6: Herwig 6 |l2l US] subsubdirectory. It will be used for the future 
development. 

— HERWIGPP: Herwig-|-- 1- [HI HS] subsubdirectory. It will be used for the 
future development. 

— interface: Some interface hies are included in this subsubdirectory. For exam¬ 
ple, QEDPS_interface.f90 is a hie to interface HELAC-Onia with QEDPS. 

• analysis. A subdirectory for perfroming analysis. 

— hbook: Hbook hies (a simplihed version written by M. Mangano) for plotting. 

— user: user dehned plot hies,like plot_user.f90. Some examples are also given. 

— PYTHIA8: the analysis code for generating histograms or Root trees by 
using Pythia 8 and Fast Jet [38] (or its core functionality FJCore). 

— heptoptagger: the HEPTopTagger [39l SOI E] source code for top quark 
tagging in the analysis stage. 

— include: some including hies for example HEPMC90.INC for dehning HepMC [M] 
common variables. 

— various: some useful tools at the analysis stage. 
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— TMVA; some examples for using TMVA contained in Root for multiply vari¬ 
able analysis. 

— LesHouches: some useful tools for dealing the Les Houches event files. 

— HepMC: a code to convert HepMC [M] file to histograms or Root trees by 
using Fast Jet or FJCore. 

• jets. A subdirecotry containing jet related tools. 

— fastjet: code for interfacing FastJet to HELAC-Onia. 

— fjcore: the source code of FJCore as well as the interface code to HELAC- 
Onia. 

— merge; the different multiplicity leading-order matrix elements and parton 
shower merging code. 

• cernlib. A subdirectory containing cernlib hies. 

— minuit: Minuit |1U6] source hies. 

• decay. A subdirectory for decaying hnal state particles. 

— decayJist.txt: a list of available decay processes. 

— Decay-interface.f90; the main decay hie. 

— DecayInfo.fOO; a hie to read the decay information from decay_user.inp in 
input subdirectory. 

— HOVll.fQO: the angular distribution hie for a vector decays into two leptons. 

— HO_chi2psia.f90: the angular distribution hie for y particle decays into a 
= 1 quarkonium and a photon. 

— HO_t2bw.f90: the angular distribution hie for top quark decays into a bottom 
quark and a W boson. 

• cluster. A subdirectory containing Python scripts. 

— create_subdir.sh; a bash shell script for creating subdirectories in the working 
directory. It is useful for running on the cluster or in the multi-core mode. 

— bin: A subsubdirectory that contains executable script hie ho_cluster after 
conhgurating and make. 

— pythoncode: A subsubdirectory that contains the python source codes. 
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* cluster.py: A file includes various cluster classes. 

* misc.py; Helpful functions defined to perform routine administrative I/O 
tasks. 

* coloringJogging.py: A file with logging color. 

* extended_cmd.py: A file conntaining different extension of the cmd basic 
python library. 

* files.py: A hie contains useful classes for dealing with hie access. 

* helaconia_run_interface.oy and helaconiaJnterface.py: A user-friendly 
command line interface to access HELAC-Onia features. 

• addon. A subdirectory for some ad hoc codes. 

— addon_process.dat: A list of available addon processes. 

— pp_psipsi_DPS: An ad hoc code for DPS of pp{p) —)■ Q 1 Q 2 + where Qi = 
J/i;,ij{2S),T{lS),T{2S),T{3S). 

— pp_psiX_CrystalBall: An ad hoc code for pp{p) —^Q + X via crystal ball 
function, where Q = J/'ip,'tp(2S),T(lS),T(2S),T(3S), Xco, Xa, Xc 2 and XbjinP) 
with J = 0,1,2,n = 1,2, 3. 

— fit_pp_psiX_CrystalBall; An ad hoc code for htting crystal ball function to 
the experimental data of pp(p) —)■ Q -|- X, where Q = J/ijj, 4’{2,S). 

— fit_pp_upsilonX_CrystalBall: An ad hoc code for htting crystal ball function 
to the experimental data of pp{p) Q + X, where Q = T(IS'), T(2S'), T(3S'). 

— fit_pp_QQ_CrystalBall: An ad hoc code for htting crystal ball function to 
the experimental data of pp{p) Q + Q, where Q is charm or bottom quark. 

— pp_QQ_CrystalBall; An ad hoc code for generating events of pp{p) Q + Q 
via crystal ball function. 

— pp_aajj_DPS: An event generator for producing pp(p) —)■ 77 -l-dijet from DPS. 

There are other subdirectories under the main directory. All generated libraries will be 
put in lib subdirectory. All module (object) hies will be put in mod (obj) subdirectory. 
Executable hies will be generated in the subdirectory bin. 


B Particle symbols in HELAC-Onia via Python script 


In this appendix, we will introduce the new particle symbols for using HELAC-Onia 2.0 
with Python scripts. We list them explicitly in Tabs,3|4|[5|6 
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Particle 

Particle ID 

Particle Symbol 

z/e, e“, M, d, /i", c, s, Ur, r~,t,b 

1,.. 

.,12 

ve,e-,u,d,vm,m-,c,s,vt,tt-,t,b 

ve~, e+, u~, d~, vm~, m+, 

Ue, e+, u, d, Uf,, /i+, c, s, ]7r, i, b 

-1,.. 

.,-12 


^,Z,W+,W-,g 

31,. 

..,35 

a,z,w+,w-,g 

H, X, 0^, <t>~ 

41,. 

..,44 

h,g0,g+,g- 


Table 3: The identity numbers and symbols of the SM “elementary” particles in HELAC- 
Onia 2.0. 


Fock State Particle ID Particle Symbol 



441001 

cc~(lS01) 


441008 

cc~(lS08) 


443011 

cc~(3Sll) 

cc[^4®'] 

443018 

cc~(3S18) 

cc[^p|^'] 

441111 

cc~(3Pll) 

cc[Vf’] 

441118 

cc~(3P18) 

CcW=0,l,2] 

4431J1 

cc~(3PJl) 

ccfPtox^] 

4431J8 

cc~(3PJ8) 


Table 4: The identity numbers and symbols for the charmonia in various Fock states in 
HELAC-Onia 2.0. 
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Fock State 

Particle ID 

Particle Symbol 

66^4''] 

551001 

551008 

553011 

553018 

551111 

551118 

5531J1 

5531J8 

bb~(lS01) 

bb~(lS08) 

bb~(3Sll) 

bb~(3S18) 

bb~(3Pll) 

bb~(3P18) 

bb~(3PJl) 

bb~(3PJ8) 


Table 5: The identity numbers and symbols for the bottomonia in various Fock states in 
HELAC-Onia 2.0. 


Fock State 

Particle ID 

Particle Symbol 


451001 

cb~(lS01) 

c6[‘sf] 

451008 

cb~(lS08) 


453011 

cb~(3Sll) 


453018 

cb~(3S18) 


451111 

cb~(3Pll) 

c6[^pf’] 

451118 

cb~(3P18) 


4531J1 

cb~(3PJl) 

rtlhSloaJ 

4531J8 

cb~(3PJ8) 

tel'si''] 

-451001 

bc~(lS01) 


-451008 

bc~(lS08) 


-453011 

bc~(3Sll) 


-453018 

bc~(3S18) 

6c[^pf’] 

-451111 

bc~(3Pll) 

6c[^Pf'] 

-451118 

bc~(3P18) 


-4531J1 

bc~(3PJl) 


-4531J8 

bc~(3PJ8) 


Table 6: The identity numbers and symbols for the mixed flavour quarkonium family 
in various Fock states in HELAC-Onia 2.0. 
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C New input parameters 


Some of the parameters in input/default.inp and input/user.inp have already been 
introdnced in Ref. [18]. The new parameters we introdnced in the new version are: 

1. energy.beainl and energy_beam2 are the energies in nnit of GeV of the hrst beam 
and second beam respectively. 

2. f ixtarget is a flag to compute the cross section in a hxed-target collision envrioment 
(T) or not (F). 

3. ranhel is a parameter to determine whether the program uses the Monte Carlo 

sampling over the helicity conhgurations. In HELAC-Onia 2.0, we extend ranhel 
to be 4, which is at the same level of performing Monte Carlo over the helicity 
conhguration with ranhel=3. Instead of using ^ d0e^(e^)*) to perform the helicity 
summation where = Yl\=± o we select the helicity eigenstate of external 

particle when ranhel=4 to take a subsequent spin-entangled decay. 

4. pdf is the PDF set number proposed in LHAPDF [6l| or in pdf/pdflist.txt. En¬ 
tering 0 means no PDF is convoluted. If one wants to use LHAPDF, please edit 
input/ho_configuration.txt and set the parameter lhapdf to be T. 

5. reweight.pdf is a flag to use reweighting method to get PDF uncertainty. It 
only works when using LHAPDF.Correspondingly, one should also specify the hrst 
(pdfjnin) and the last (pdfnnax) of the error PDF sets. 

6. reweight.Scale is a hag to use reweighting method to get renormalization and 
factorization scale dependence, which requires alphasrun=T. One can change the 
lower bound and upper bound for renormalization/factorization scale variations via 
parameters rw_RScale_down,rw_RScale_up,rw_FScale_down and rw_FScale_up. 

7. useMCFMrun is a hag to perform the strong coupling as renormalization group run¬ 
ning in the MCFM [101] way. 

8. toodrawer.output, gnuplot.output, root.output are hags to ask HELAC-Onia 
to plot histograms and to output into TopDrawer, Gnuplot and Root hies on the 

fly- 

9. emep_ISR_shower is a parameter to determine whether use QEDPS to take into 
account initial state radiation ehects in electron-positron collisions (1) or not (0). 
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10. parton_shower is a parameter to determine whether perform parton shower. parton_shower=0 
means no shower, i.e. hxed-order calculation. The shower can only be used when the 
corresponding parton shower program is already installed and the user has already 
edited properly in input/ho_configuration.txt. 

All other parameters are listed in default. inp. The user can £x his/her values in 
user. inp following the format in default. inp. 


D Addon codes 

In this section, we will describe some addon codes implemented in HELAC-Onia 2.0 for 
dedicated studies. All of the addon codes have been listed in addon/addon_process.dat. 


D.l Single-quarkonium hadroproduction with crystal ball func¬ 
tion 


In fact, the description of the quarkonium-production mechanisms is still a challenge for 
theorists, especially current the state-of-the-art computation in NRQCD cannot describe 
the single-quarkonium-hadroproduction data in the whole kinematical region. It would 
be quite interesting and might be necessary to use emiprical function to describe the 
single-quarkonium production from pp or pp collisions with a data-driven way and use 
it to test other mechanisms like DPS or pA and AA collisions. Moreover, it also pro¬ 
vides an economy way to generate events for single-quarkonium hadroproduction. There¬ 
fore, HELAC-Onia 2.0 has already been implemented some dedicated codes to £t the 
single-quarkonium-hadroproduction data and to generate events of the single-quarkonium 
hadroproduction. 

Let us start with the description of the calculation. The initial-averaged squared am¬ 
plitude for single-quarkonium Q hadroproduction with the assumption of the dominance 
of gluon-gluon channel can be expressed in a crystal ball function [H] 


\Agg^Q+x\‘^ = 


Kexp{-K^] 


Wlltill Jrj' ^ \-i T/ 


+ when Pt>{Pt} 


( 16 ) 


where K = X^ks/Mq and s is the partonic center-of-mass energy. Then the cross section 
of single quarkonium Q production in pp collisions is 


a{pp Q + X) = J dxidx2fg{xi)fg{x2)^\Agg^Q+x\‘^dLIPS, (17) 
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where fg is the gluon PDF and dLIPS is the Lorentz-invariant phase space measure for 
PP — t Q + X. The coefficients A,K,n and (Ft) can be determined by htting it to the 
experimental data. 

D.1.1 Fit codes 

The codes for fitting '0(15', 25') and T(15', 25', 35') are in subdirectories fit_pp_psiX_CrystalBall 
and fit_pp_upsilonX_CrystalBall respectively. One can use the command lines 

I H0> generate addon 3 

and 

I H0> generate addon 4 

to drive the corresponding codes to perform htting with Minuit |lU6j package. The input 
hies dedicated to these ad hoc codes are in fit_pp_psiX_CrystalBall/input and 
fit_pp_upsilonX_CrystalBall/input. One can specify the meson in state.inp and the 
htting parameters in fit_param_card.inp. The selected experimental data can be as¬ 
signed in dataJistJ.inp, where i is the number in state.inp. Some htted results are con¬ 
tained in fit_pp_psiX_CrystalBall/fitresults and fit_pp_upsilonX_CrystalBall/fitresults. 
We have checked the htted results of Ref. [9] for prompt J/il) production with the same 
setup. 

For instance, through a combined ht of d'^a/dPTdy to the ATLAS [70], CMS mm, 

LHCb |lU8j and CDF |1U9] data, we obtained k = 0.543 and A = 0.118 for prompt '0(25') 
when (Ft) = 4.5 GeV and n = 2, where = 242 for total 90 experimental data. The 
comparisons are shown in Figj^ The result is collected in fitresults/psi2s/fitl. 

D.1.2 A simple event generator for single-quarkonium hadroproduction 

With the htted parameters, we wrote a simple event generator for pp{p) Q + X, where 
Q = J/0,0(25'),T(15'),T(25'),T(35'),yco,Xci,Xc2 and Xbj{nF) with J = 0,1,2,n = 

1,2,3. The code is put in the subdirectory addon/pp_psiX_CrystalBall. One can 
drive such program with the following command line 

I H0> generate addon 2 

Some special input parameters can be specihed in pp_psiX_CrystalBall/input. One 
can set the type of Q in state.inp and its polarization in polarization.inp. The hie 
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crystalball.inp is used to input the parameters X,K,n and (Pt)- We have performed 
some applications in Ref. m; 


D.2 Double parton scattering for double-quarkonium produc¬ 
tion 


We also implemented the code for calculating DPS for double-quarkonium production in 
a pp or pp collider. One of its application can be seen in Refs. [iniiiii]. We used a simple 
but widely-used “pocket formula” to describe DPS for double-quarkonium production 
PP Q1Q2 + w 


a 


DPS _ 
S1S2 


1 + ^QiQ 2 ^eff 


(18) 


where ctq. is the cross section for single quarkonium Qi production and aes is a parameter 
to characterise an effective spatial area of the parton-parton interactions. cTefr can be 
related to the parton spatial density F{h) inside the proton as 




j d"b(F(b))" 


(19) 


Within the factorization, aes should be independent of hnal states but it might change 
with different species of initial partons and its Bjorken fraction x. A hrst order assump¬ 
tions of (Jeff is independent of process and energy, which however should be checked case 
by case. 

For single-quarkonium production, we used the crystal ball function described in ap¬ 
pendix |D.1| to estimate the squared amplitude. This special code can be found in ad- 
don/pp_psipsi_DPS. The command line to generate this process is 

I H0> generate addon 1 


Similar to other addon codes, the input parameters dedicated to this code is in pp_psipsi_DPS/input. 
The parameter Ueff can be specified in sigma_eff.inp. The user can change the type of the 
quarkonium pair Qi and Q 2 in the hie states.inp. The polarizations of Qi can be hxed 
in polarization_<name>.inp, where <name> is the name of Qi,i.e. J/tl) = jpsi, ■^(25') = 
psi2S,T(lS') = Y1S,T(2S') = Y2S,T(3S') = Y3S. The parameters A,K,n and {Pt) in 
the crystal ball function should be told in the hies crystalball_<name>.inp. 
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D.3 Heavy-flavor quark pair hadroproduction with crystal ball 
function 

Using the crystal ball function for heavy-flavor quark pair or open heavy-flavor meson pair 
production, one can also perform a £t to experimental data at the hadronic colliders in 
the Pt spectrum of the quark/meson production. Hence, it also provides an opportunity 
to use a data-driven way to analyze the corresponding heavy-flavor quark/meson pair 
production, which usually suffers large theoretical uncertaities in a perturbative com¬ 
putation. Such a method indeed has been applied in the open charm production at a 
proposed hxed-target experiment at the LHC (AFTER@LHC) in Ref. [HD]. The £t can 
be performed using the following commands 

I H0> generate addon 5 

Some input parameters for fitting are needed to be specihed in fit_pp_QQ_CrystalBall/input. 
Moreover, with the htted parameters, one can use the 

I H0> generate addon 6 

to generate the unweighted events for the heavy quark pair production in proton-proton 
or proton-antiproton collisions. 


D.4 


Double parton scattering for associated production of dipo- 
ton and dijet 


Similar to Eq.(18), we have a “pocket” formula for diphoton and dijet production via 
DPS mechanism in pp collisions 


DPS 

‘^77+ii 


O’-yyO'jj ^ <^7+j‘^7+j 

^efF 2(Tgff 


( 20 ) 


The LO matrix elements of a 6 —)■ 7 -|- 7,7 -f j,j -|- j have been implemented in HELAC- 
Onia with the correct color flow. To the diphoton production, we also implemented 
the gluon-gluon initial state process, which is a loop-induced process for diphoton pro¬ 
duction. However, due to the high luminosity of the gluon-gluon initial state at a high 
energy collider, such a contribution might be substantial. Unweighted events for the DPS 
contributions to diphoton and dijet production can be generated by the command 


I H0> generate addon 7 
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One can change the input aes in the input file pp_aajj_DPS/input/sigma_eff.inp. 
There is also a file pp_aajj_DPS/input/subprocess.inp for user to select/drop some 
partonic subprocesses and to choose to generate the unweighted events of pp —)■ 7 + 7 , pp —)■ 
7 + j,pp —)■ j + j instead of the DPS process. Such a functionality is very useful to cross 
check and to specify a global K-factor from the missing higher-order quantum corrections. 
Some studies on the inclusive DPS production rates of this process at the Tevatron |112j 
and the LHC |113j have been explored in the literature. 
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Figure 2: The differential distributions for pp —)• J/^|J + J/^|J + cc in the CMS fidicuial 
region [73]: (a) absolute azimutal difference; (b) absolute rapidity difference ; (c) invariant 
mass distribution; (d) the vectorial transv^e momentum sum; (e) leading pp] (f) sub¬ 
leading pt- 














































































Figure 3: Illustrative plots for J/'ijj production at 13 TeV LHC with parton shower from 
Pythia8.186. We presented the hxed-order LO calculation (solid curve), LO+PS (long- 
dashed curve), LO-I-PS but turning off ISR (short-dashed curve) and LO-I-PS but turning 
off FSR (dotted curve). 
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Figure 4: Validation of lepton angular distributions in —)■ . 
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Figure 5: Validation of J/'0 angular distributions in Xci J/ijj -f 7 . 
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Figure 6: Validation of J/V’ angular distributions in Xc 2 J/V’ + 7- 
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Figure 7: Program structure of HELAC-Onia with version 2.0. 
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Figure 8: Combined fit of cPa/dPTcly to ATLAS [TO] (lst-3rd plots), CDF |1U9] (4th 
plot), CMS |lU7j (5th-7th plots), LHCb jlU8j (8th plot) for prompt '0(2*5') production. 
The plots are generated automatically by HELAC-Onia 2.0. 
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